/*  ldbl2ary.cpp  by K.Tsuru  */
/*****************************
SN class's auxiliary function
long double --> DRADIX conversion
******************************/
#ifndef SN_H
#include "sn.h"
#endif

/*
It gets the values of "e" and "d" which satisfy a condition
p = d * 10^e (0.1 =< d < 1.0 ).
Note that even if p = 1.0*10^e, *d = 0.999999999... and becomes one by adding
DBL_EPSILON to that.
*/
static int dblExp(ldouble p, ldouble* d){
    ldouble absp = fabs(p); // check |p| == 1.0 or 0.1
    if(absp == 1.0){        // 1.0 = 0.1 * 10^1
        *d = p > 0.0 ? 0.1 : -0.1;  return 1;
    } else if(absp == 0.1){ // 0.1
        *d = p > 0.0 ? 0.1 : -0.1;  return 0;
    }

    ldouble q = log10(absp);
    int e = (int)q, e_add = 0;
    if(q >= 0.0) e++;
    if(e > DBL_MAX_10_EXP){  absp /= 10.0;  e--;    e_add = 1;  }
    if(e) q = absp*pow(10.0, -e);
    else  q = absp;

    if( q + DBL_EPSILON >= 1.0){// q = 0.99999999...*10^e = 0.1 * 10^(e+1)
        q = 0.1;    e++;
    } else if( q < 0.1){    // q = 0.09999999.... *10^e = 0.1*10^e
        q = 0.1;
    }
    *d = p > 0.0 ? q : -q;
    return e + e_add;
}
/*
It converts "d" to fType array "ary" where 0.1 =< d < 1.0 and return the size.
*/
static uint ConvertfType(FigBlock& ary, ldouble d){
    int j, len = 0, r;
    fType under = 0, round = 0;
    for(j = 0; (d > 0.0) && (len <= DOUBLE_FIG); j++){
        ary[j] = (fType)d;
        d -= (ldouble)ary[j];
        d *= (ldouble)DRADIX;
        len += DFIGURES;
    }
    ary[j] = (d >= 1.0) ? (fType)d : 0;
    r = len - DOUBLE_FIG;
    if(r > 0 && ary[j]){
        under = 1;
        while(r){   // r > 0 then under >= 10
            under *= 10;    r--;
        }
        //Check last "r" figures.
        round = ary[j] % under;     // ary[j] = 8999, under = 10, round = 9
        if(round >= under/2 ){  // 9 > 10/2 = 5
            ary[j] -= round;    // ary[j] = 8990
            ary[j] += under;    // ary[j] = 8990 + 10 = 9900
        } else round = 0;
    }
    if(round){
        int k;
        for(k = j; k > 0 ; k--){    // normalize
            if(ary[k] < DRADIX) break;
            else {
                ary[k-1]++;
                ary[k] = 0; //len -= DFIGURES;
            }
        }
        if(k < j) under = 0;
        j = k;
    }
    if(under){  //cut out last figure
        ary[j] /= under;    ary[j] *= under;
    }
    while(!ary(j) && j) j--;//some low figures becomes zero by round
    if(ary(0)) return 1;    //d = 0.9999999999 -> d = 1.0 ( ary(0) = 1 )
    return (uint)j+1;
}
/**********************************************************
main body
Provides radix conversion ldouble "x" to res[] in DRADIX.
It returns the value of exponent.
***********************************************************/
int ldoubleToArray(ldouble x, FigBlock& res, uint *size){
    res.reserve(minArraySize-1);
    res[0] = res[1] = 0;
    *size = 2;
// x == 0.0 ?
    if(x == 0.0) return 0;

    ldouble absX = fabs(x);
    int exp;
    ldouble     ip, dp;
    dp = modf(absX, &ip);
    if(dp + DBL_EPSILON >= 1.0){
        dp = 0.0;   ip = ip + 1.0;
    }
//An integer which has an order of "DRADIX" can be exactly represented by ldouble.
//Then "+ROUND_DBL_INT" is not nesessary.
    if( (dp == 0.0) && (ip < (ldouble)DRADIX ) ){//integer of one figure
        res[1] = (fType)ip;
        return 1;
    }

    res.reserve(DOUBLE_FIG/DFIGURES+3);
    ldouble d;
    exp = dblExp(absX, &d); // absX = d*10^exp, 0.1 <= d < 1.0

    fType p10;
    int dexp, r;
    uint sz = ConvertfType(res, d);
    if(sz == 1){    // d = 0.9999999999 -> d = 1.0
        res[0] = 0;     res[1] = DRADIX/10;
        sz = 2;         exp++;  //do not return to settle the value of "exp"
    }

    dexp = exp/DFIGURES;
    r = exp - dexp*DFIGURES;
    if(r > 0){  r -= DFIGURES;     dexp++;  }   //r <= 0
    p10 = 1;
    while(r < 0){ p10 *= 10;    r++;    }

    uint j;
    if(p10 > 1){    //divide by p10
        res.reserve(sz);
        res[sz] = 0;    sz++;
        ulong w, u;
        u = 0;
        for(j = 1; j < sz; j++){
            w = u*DRADIX +(ulong)res[j];
            res[j] = fType(w/p10);
            u = w%p10;
        }
    }
    j = sz-1;
#ifndef NDEBUG
res(j);
#endif
    while(res[j] == 0) j--;
    *size = j+1;
    return dexp;
}